Charge-density- wave instabilities driven by multiple umklapp scattering 
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We show that the concept of umklapp-scattering driven instabilities in one-dimensional systems 
can be generalized to arbitrary multiple umklapp-scattering processes at commensurate fillings given 
that the system has sufficiently longer range interactions. To this end we study the fundamental 
model system, namely interacting spinless fermions on a one-dimensional lattice, via a density matrix 
renormalization group approach. The instabilities are investigated via a new method allowing to 
calculate the ground-state charge stiffness numerically exactly. The method can be used to determine 
other ground state susceptibilities in general. 



I. INTRODUCTION 

The relevance of umklapp-scattering terms in one- 
dimensional systems with commensurate filling has long 
been established,* Besides the general theoretical interest 
in one-dimensional systems because of the existing ex- 
act solutions to some of thern^ and the field-theoretical 
description of the low energy physics^ the application 
of one-dimensional models to the rich phase diagrams 
of quasi one-dimensional systems such as the Bechgaard 
and Fabre salts^ has a long standing history. The de- 
generacy of the ground state leads to the presence of 
solitonic excitations in the one-dimensional systems and 
recent studies focus on the question to what extent these 
collective modes are observable 5 - in the systems under in- 
vestigation. 

Here we focus on the generic instabilities of systems 
of fermions without the internal magnetic degree of free- 
dom. The Hamiltonian under consideration is 
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with fermionic creation (cj) and annihilation (c, ) op- 
erators at site I, nearest- neighbor hopping amplitude 
t, density-density interaction amplitude V s between 
fermions on sites of separation s and density operators 
ni = cjc l . Energies are given in units of t. The de- 
cay of the interaction with distance requires in general 
V s +i < V s , with the exception of zig-zag chains and 
some spin systems where the nearest-neighbor interac- 
tion is suppressed by the lattice geometry. We apply pe- 
riodic boundary conditions, i.e., cm+i = en* The model 
in Eq. JIJ is directly applicable to spin 1/2 chains that 
can be mapped onto chains of spinless fermions^ Fur- 
thermore, its widely studied continuum representation 8 
is similar to the charge sector of the bosonized Hubbard 
model. 5-9 Finally, Eq. allows to study the interplay of 
ordered phases with different modulation wave vectors, 
which proves useful for the understanding of materials 
that exhibit multiple phase transitions^ 

The purpose of this paper is threefold, (i) We present a 
novel method to determine the ground-state curvature or, 
equivalently, the ground-state charge stiffness, (ii) The 



approach allows for an accurate identification of charge- 
density-wave (CDW) instabilities at commensurate fill- 
ings. The phase diagrams for various relevant sets of pa- 
rameters are shown, (iii) We discuss the physical insight 
that can be gained from the continuum representation 
of Eq. JQ) and by comparison with the numerical results 
show the importance of lattice effects. 

To this end we present a detailed density matrix renor- 
malization group (DMRG) studyii of the lattice model 
Eq. QJ. Sufficiently large system sizes are accessible to 
model the thermodynamic limit without the shortfall of 
neglecting possible lattice effects. 



II. NUMERICAL APPROACH 

The ground state curvature C is defined as 
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C = M 
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where is a magnetic flux that penetrates our system 
which is closed to a ring. The magnetic flux inside the 
ring can be gauged into a single bond leading to twisted 
boundary conditions c\ — e^CM+i- According to Kohn* 2 
C is equivalent to the Drude weight of the T = dc 
conductivity. The curvature in Eq. J2J is the second order 
coefficient of the Taylor series of Eo(<j>), i.e., Eq(4>) = 
E o (0) +EW<f>+ 2^(/> 2 +0((f> 3 ). In order to compute the 
coefficient exactly we expand H(4>) and the ground-state 
wave function \^q((/))) at ^ = in second order. H is 
given via 



^)=i/(°)+z0J-^-T + O(0 3 ), 
where = iJ is the current operator with 
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and 



T = H^ = -{c\ [Cl +c\c M ) 
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is the kinetic energy on the bond between the Mth and 
first site. Due to the normalization of the wave function 
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(*(</>) | = 1 = | the first order correc- 

tion I^W) in the expansion | *((/>)) = |* (0) ) + </>|^ (1) ) + 
0(4> 2 ) is orthogonal to the ground state l^ )) at zero 
flux and is obtained by solving the set of linear equations 

(HW - 4 0) ) |*«) = (V> - tj) |*<°>) . (6) 

on the subspace that is orthogonal to the ground state. 
The resulting curvature is given by 

C = M(¥V\T\¥V) + 2M(^\iJ - E^\¥^) (7) 

with £?w = (*(°)|irW|*( )). 

It must be stressed that this approach is exact. There- 
fore the accuracy of our numerical results is only limited 
by the truncation error of the DMRG. Consequently any 
ambiguity in the determination of C incurred through a 
numerical finite difference approximation of the deriva- 
tive of Eq{4>) is avoided. Moreover, our procedure allows 
to calculate C without referring to a system with a finite 
flux that could destroy the commensurability effects we 
want to measure. 

For systems with time-reversal symmetry and a non- 
degenerate ground state one finds i?W = and from Eq. 
© follows that l 1 !^ 1 )) is purely imaginary. In turn it fol- 
lows that C can be calculated by using real numbers only. 
Therefore solving the linear set of equations © and stor- 
ing the additional target vectors ImlvE^ 1 ))) and J\^ 1 ')) 
are outweighed by the fact that we can restrict ourselves 
to real numbers. Since 

(^f(i) | $(0)) = 0j solving Eq. © 
is stable and can be performed by a mixture of standard 
preconditioned minimal residue and conjugate gradient 
methods which we extended by a projection on the space 
orthogonal to the subspace of the ground state(s). The 
extension to degenerate ground states is straightforward 
and can be performed similar to Brillouin-Wigner per- 
turbation theoryiH For the small systems (M < 30) we 
keep 300 to 400 states per block A and B, while we keep 
up to 1000 states per block for the larger systems. In ad- 
dition to the ground state l^' ') and the auxiliary states 
J\^^) and l^ 1 )) we target at least for the first four 
excited eigenstates to treat degeneracies properly. 

III. PHASE DIAGRAMS 

In order to investigate the general features the model 
in Eq. (JJJ we study its phase diagram as a function of 
the nearest- (Vi) and next-nearest-neighbor (V2) interac- 
tion. The result is presented in Fig. ^ where the curva- 
ture is plotted versus V\ and V2. Charge order implies 
(7(Vi, V2) — 0. In Fig. IHa) the system is shown at half 
filling for a chain of length 32. In CDW phase I the 
ground state is twofold degenerate^ with ordering pat- 
tern (»o»o) and (o«o«). Here o denotes a vacant and • 
denotes an occupied site. In phase II the ground state is 
fourfold degenerate with ordering pattern (••00), (o»»o), 
(oo»»), and («oo«). The transition between phases I and 




FIG. 1: Curvature of the (a) half and (b) one-third-filled 
chain with (a) M = 32 and (b) M = 24 sites as a function 
of the nearest- (Vi) and next-neares-neighbor (V2) interaction 
strength. Charge order implies (7 = 0. The ordering pattern 
are discussed in the text. In (a) the CDW phase I doubles 
the unit cell, phase II quadruples the unit cell and the ground 
states are twofold and fourfold degenerate, respectively. In (b) 
the ordered phase triples the unit cell with triply degenerate 
ground state. 



II along the line V± — 2V2 is expected from simple poten- 
tial energy arguments for the corresponding real-space 
ground-state configurations. The finite width of the uni- 
form region is due to a gain in kinetic energy near the 
line V\ = 2V2, where the charge ordering effects are re- 
duced by the competition of the two potential energies. 
Figure^b) shows the curvature of a chain of length 24 at 
1/3 filling. The ordered phase is observed for sufficiently 
large values of V\ , V2 > 4 and has a triply degenerate 
ground state with ordering pattern (•oo»oo), (o»oo»o), 
and (ootoot). 

Figure[21shows the finite size effects of the curvature at 
half filling via cuts through Fig. ^a) for various system 
sizes. FigureEJa) shows C(Vi,0), panel (b) C(0, V2), and 
panel (c) C{V U 10 - 2V*i). For V 2 = and M — > 00 it 
is known from Bethe anstatz calculations^ that in the 
ordered phase for V\ > Vi, c = 2 an excitation gap opens 

exponentially 16 as ~ cxp ^— 7r 2 / ^/8(Vi — Vi jC )^ . As a 

result the convergence of the curvature to the thermody- 
namic limit in CDW I near the transition point is very 
slow [Fig. |2Ja) and Ref. 17]. In contrast we find a gap in 
CDW II obeying the power law A n | - (V 2 - V 2 , c ) 3/4 
as shown in the inset of Fig. Efc). The reason for 
the different behavior lies in the bond-ordered phase 
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FIG. 2: Curvature at half filling as a function of the interac- 
tion parameters for different system sizes along cuts through 
Fig. 13 a). Panel (a) V 2 = 0. Panel (b) Vi = 5 - 2Vi. Panel 
(c) Vi = 0, the inset shows that the gap in CDW II opens as 
An ~ (V2 — V2,c) 3/ ^ 4 - Panel (d) shows a sketch of the resulting 
phase diagram. The dashed line indicates the direction of the 
cut shown in panel (b) , the dash-dotted lines shows 2 V2 = Vi . 
Dots mark the transition points determined herein. 
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FIG. 3: Level crossings at the phase boundaries between the 
BO and CDW II phase and changes of sign of the finite size 
effects of M{E\ — Eo) at the boundaries of the LL phase along 
the line V2 = 5 — 2Vi. System sizes are M = 24 (diamonds), 
M = 36 (crosses) and M = 48 (triangles). Lines are guides 
to the eye. 



numerically within this work. The numerical analysis for 
large values of V* = 20, 50, 100 with V 2 = V* - 2V X sug- 
gests that for V* t the uniform region (BO and LL 
phases) between CDW I and CDW II extends along the 
line V\ = 2V2 with a width of ~ \/5£- 19 At sufficiently 
high interaction strength the LL phase disappears. The 
triple point between the BO, LL, and CDW I phase can 
be estimated to be near (Vl,t,V2,t) ~ (8.2,3.6), see also 
Sec. El 



(BO) equivalent to that found in frustrated Heisenberg 
chains, 18 which we find strong evidence for between the 
Luttinger liquid (LL) and CDW II phases, see Fig.[2Jd). 

The numerical determination of the critical values 
of the interaction is difficult because of the smallncss 
of the gap near the transitions. It turns out though 
that the finite size effects of M{E\ — Eo) change sign 
at the phase transition between the LL and the or- 
dered BO and CDW I phases, i.e., liniM^oo dM(E\ — 
E )/(dM)\v 1 =v ltC = 0. On the other hand, in the 
CDW II phase systems with sizes Afmod4 = have much 
smaller finite size effects than for Mmod4 = 2 while in 
the BO phase both are equivalent. Moreover, the systems 
exhibit well defined level crossings between the bond or- 
dered and the CDW phases^ The effects are illustrated 
in Fig. |3 for M = 24 (diamonds), M = 36 (crosses), and 
M = 48 (triangles) for a cut along the line V2 = 5 — 2Vi 
for the four lowest excited states. Lines are guides to the 
eye. Note that for the actual parameter determination 
systems of sizes up to 60 sites where studied near the 
phase transitions. 

Combining different approaches we estimate for V\ = 
a value of V2, c = 2.66 ± 0.1 while along the line V2 — 
5 — 2Vi the critical points into the CDW phases are 
(Vi,c,V 2 ,c) € {(1-05 ± 0.05,2.9 ± 0.1); (2.4 ± 0.05,0.2 ± 
0.1)}. The transitions from the LL to the BO phase occur 
at (Vi, b , V 2 , b ) = (0,2.15±0.1) and (1.35±0.05, 2.3±0.1). 
The resulting phase diagram is sketched schematically in 
Fig. |2Td) . Thick dots are transition values determined 



IV. COMMENSURABILITY EFFECTS 

Any deviation from commensurate filling in one- 
dimensional systems leads to the presence of solitonic 
excitations even at zero temperature which disorder the 
system due to the degeneracy of the ground stated Our 
approach allows to demonstrate the sensitivity of our one- 
dimensional model system to the commensurability of the 
filling very nicely as is shown in Fig. The diamonds 
show the curvature for different fillings at different sys- 
tem sizes 18 < M < 60 for V x = 10 and V 2 = 3, i.e., the 
system is charge ordered both for fillings of 1/2 and 1/3 
as can be seen from Figs.^a) and^b). Since the values 
of interaction are such that the system at N/M = 1/3 
is close to the transition to the non-ordered phase the 
convergence is slow at that point. The squares show the 
curvature for V± = 10 and V% = 5, where no charge or- 
dering is observed at half filling due to the competition 
between Vi and V 2 - Note that the results are symmetric 
about N/M — 1/2 because of particle-hole symmetry. 

The results for Vi, V2 7^ discussed in detail above can 
readily be generalized to systems with longer range in- 
teraction. For example, for Vi, V%, V3 ^0 and sufficiently 
large V3 we find a phase with ordering pattern (»ooo) and 
four-fold-degenerate ground state at quarter filling. Since 
a detailed discussion would be a simple extension of the 
case studied above it is omitted here. The crucial result 
that follows from the generalization is that for all sys- 
tems under investigation (s max = 1, 2, 3, 4, 5) an instabil- 
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FIG. 4: Curvature as a function of filling for different num- 
bers of fermions 2 < N < M/2 at system sizes 18 < M < 60. 
Diamonds: V\ = 10 and V2 = 3 with charge order at 
N/M = 1/2 and N/M = 1/3. Squares: Vi = 10 and V 2 = 5 
with charge order only at N/M = 1/3 (c.f. Fig. 0. 




FIG. 5: Curvature for s max = 5 as a function of filling for 
different numbers of fermions 2 < TV < M/2 at system sizes 
20 < M < 48 with Vi = 120, V 2 = 60, V3 = 30, V 4 = 15, 
and V5 = 6 with charge order instabilities at N/M — 1/6, 
7V/M = 1/5, N/M = 1/4, 7V/M = 1/3, and N/M = 1/2. 
Note that for the strong interactions under investigation there 
is also a second-harmonic anomaly at N/M = 2/5. The result 
supports the notion that instabilities occur only at fillings 
with N/M > l/( 



ity at commensurate filling N/M is only observed for suf- 
ficiently long-ranged interaction, i.e., V s ^ V s < s max 
with s max = M/N—l. The case of s max = 5 is illustrated 
in Fig. for V x = 120, V 2 = 60, V 3 = 30, V4 = 15, and 
V5 = 6 and exhibits charge density instabilities only for 
N/M > 1/6. Note that for sufficiently strong interaction 
there are also second-harmonic anomalies such as the one 
at N/M = 2/5 in Fig.0 

The notion of instabilities only for sufficiently long 
ranged interaction is quite intuitive when considering the 
system under investigation in real space, because the sys- 
tem cannot order when the mean distance between par- 
ticles is larger than the interaction range. 



V. CONTINUUM REPRESENTATION 

The low energy properties of the Hamiltonian Eq. 
|JT)| can be approximated by its field-theoretical equiv- 
alent obtained by the standard Abelian bosonization 
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FIG. 6: Numerical estimate nC/v for the LL parameter K for 
cuts along V2 = 5—2Vi (circles), V% = 10— 2V\ (triangles), and 
V2 = 20 - 2Vi (diamonds). For 0.5 < tyC/v = K the system 
is in the LL regime, for 0.5 > nC/v in an ordered phase. 
For Vi ~ 8 the system is close to the tricritical point where 
the LL phase disappears for larger interaction strength. The 
numerical precision of the method to determine irC/v can be 
estimated at V\ — 2 and V2 = and amounts to about ±5%. 



tcchniqueA^ The kinetic and forward-scattering parts 
lead to the free Hamiltonian 



TLa 



v 

167T 



dx K IT 



1 

K 



(8) 



with excitation velocity v, Luttinger liquid parameter K 7 
Bose fields 4>{x), and their conjugate momenta Il(a;). The 
continuum coordinate is x = lim a ^o l a with the lattice 
constant a and M — ► 00. The umklapp-scattering part is 



Ti u = lim > • 



V s 



a/20 - (4fc F - G)x 



, , dx cos 

with Fermi wave vector = np/a where p = N/M < 1 
is the number of fermions per site. G is a reciprocal 
lattice vector. For sufficiently small K < 1/2 and if 
4k-p = 2n/a, Eq. (0 yields a relevant contribution to the 
total Hamiltonian TL = H.q + Ti u which then represents a 
sine-Gordon models^ Indeed, using the finite-size scaling 
of the ground-state energy per site in the LL phas o 18 i 20 
= e(oo) — where v is the velocity of the el- 

ementary excitations, and the relation C = —K we find 
K > 1/2 in the LL phase, K 1/2 at the LL-phase 
boundary, and K < 1/2 elsewhere within the numerical 
precision of the methodic of about ±5%. The latter can 
be estimated at V\ = 2 and V2 = 0. This is illustrated in 
Fig. ©for cuts along V 2 = 5-2Vi (circles), V 2 = 10-2V4 
(triangles), and V2 = 20 — 2Vi (diamonds) and confirms 
the applicability of Eqs. © and © for p = 1/2. 

The obvious shortfall of Eq. © is that the operator in 
the integrand does not depend on the range of the inter- 
action s as a direct consequence of the continuum limit. 
The numerical study herein proves that longer range in- 
teractions lead to instabilities also for commensurate fill- 
ings other than p — 1/2 or, equivalently, for fcp ^ G/4, 
if the interaction range is larger than or equal to the 
mean distance between the fermions on the lattice, i.e., 
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for (s max + 1) > 1/p. In order to reproduce this result 
within the continuum model, Eq. @ must be modified 
to 



H u = lim 



a^o ^ 2(ira) 2 

r—l v ' 



dx cos 



\f2n r <j) — k c gx 



(10) 



with fc e ff = 2(r + l)fcp — mG. The presence of terms of 
the form of Eq. (|TU|l with A r ^ V r follows from the 
perturbative inclusion of lattice corrections which yields 
K r = (r + l)/2i2i At any commensurate filling, where 
2(r + l)kp = 2i:m/a, the term for r = m/ p— 1 is relevant 
if K < 0.5k" 2 ^ 

Preliminary numerical results suggest that K indeed 
is decreased with increasing interaction range. For ex- 
ample, near the phase transition in the case of 1/3 fill- 
ing, where s max = 2, we find K ~ 2/9 indicating that 
the term with K2 = 3/2 indeed becomes relevant as 
expected. This leads to the conjecture that the break- 
down of the LL phase is determined by a critical value 
of K = 2/ (s max + l) 2 . Further support for this result re- 
quires a detailed analysis of the dependence of the value 
of K at the LL phase boundaries as a function of s max 
which is numerically involved and beyond the scope of 
the present study. 

The excitations described by the interaction term TL U 
are r + 1 electrons that are scattered from the left Fermi 
point to the right Fermi point and vice versa. It can 
be referred to as multiple (r + 1) umklapp scattering. 
At incommensurate fillings the integrand in Eq. IjlUfl is 
oscillatory, the umklapp term is effectively averaged out, 
and no instability is observed. 



VI. CONCLUSIONS 

We have shown that (i) the exact zero-flux determina- 
tion of the ground-state curvature (or charge stiffness) is 
a powerful tool to determine numerically instabilities in 
interacting one-dimensional Fermion systems. Note that 
the method can easily be adapted to determine ground 
state susceptibilities in general, (ii) The application to 
systems of spinless fcrmions with sufficiently long-range 
density-density interaction V s , i.e., for s max > M/N — 1, 
leads to a comprehensive understanding of the rich phase 
diagrams at commensurate fillings. The approach im- 
pressively demonstrates the breakdown of long-range or- 
der for any incommensurate filling, (iii) Lattice effects 
need to be included properly in the continuum model 
in order to reproduce the numerically observed instabili- 
ties. Within the adapted field-theoretical approach they 
find an interpretation as driven by multiple (mM/N) 
umklapp-scattering processes. 
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